% This script calculates PNJunction of parabolic bands with MoS2
% parameters. Various Dopant Concentrations are considered.

% Physical Constant Setup (In SI Unit)
a0 = 0.5291772109217;  % Bohr Radius in Angstrom 
eps = 1;               % Relative permittivity
eps_0 = 8.8541878176E-12;   % Vacuum permittivity 
e = 1.6E-19; % Electronic Charge
hbar = 1.05457172647E-34;
Vh = 27.2113850560;
kB = 8.6173324E-5/27.21138386;
vf = 1E6/(2.187691263373E6);

% Energy Grid
Ne = 1000; nec = (Ne)/2;
Emax = 0.05;
Emin = -0.05;
he = (Emax - Emin)/(Ne-1);
E = linspace(Emin, Emax, Ne);

% DOS1 (In Hartree Atomic Unit)
mh1 = 0.58;
me1 = 0.48;
Eg1 = 1.85/Vh; 
Ev1 = -Eg1/2; Ec1 = Eg1/2;
D1 = zeros(1,Ne);
D1(1:round((Ev1-Emin)/he)) = 6*mh1/pi;
D1(round((Ec1-Emin)/he):end) = 6*me1/pi;

% DOS2 (In Hartree Atomic Unit)
mh2 = 0.58;
me2 = 0.48;
Eg2 = 1.85/Vh;
Ev2 = -Eg2/2; Ec2 = Eg2/2;
D2 = zeros(1,Ne);
D2(1:round((Ev2-Emin)/he)) = 6*mh2/pi;
D2(round((Ec2-Emin)/he):end) = 6*me2/pi;

% Spatial Grid (In Hartree Atomic Unit)
L = 1000;
Nx = 1000;
ncx = (Nx+1)/2;
x = linspace(-L/2,L/2,Nx);
f1 = -1E-3;
f2 = 1E-3;


%% First Step
for i = 1:1
    Q1 = f1 + sum(D1(1:round((Ev1-Emin)/he)))*he;
    Q2 = f2 + sum(D2(1:round((Ev2-Emin)/he)))*he;
    [Va, Ra, Err] = PNJunction(Q1, Q2, [E;D1], [E;D2], 300*kB, 1E-4, Nx, L);
    filename = sprintf('MoS2(f1=%1.1E,f2=%1.1E,L=%i,N=%i)',f1,f2,L,Nx);
    dlmwrite([filename, 'V.txt'],[x',Va']);
    dlmwrite([filename, 'R.txt'],[x',Ra']);
    dlmwrite([filename, 'Err.txt'], Err);
end

%% Second Step
for i = 1:1
    filename = sprintf('MoS2(f1=%1.1E,f2=%1.1E,L=%i,N=%i)',f1,f2,L,Nx);
    V0 = dlmread([filename, 'V.txt']);
    Vx = V0(:,2)';
    Q1 = f1 + sum(D1(1:round((Ev1-Emin)/he)))*he;
    Q2 = f2 + sum(D2(1:round((Ev2-Emin)/he)))*he;
    [Va, Ra, Err] = PNJunction(Q1, Q2, [E;D1], [E;D2], 300*kB, 1E-3, Nx, L, Vx);
    
    filename = sprintf('MoS2(f1=%1.1E,f2=%1.1E,L=%i,N=%i)',f1,f2,L,Nx);
    dlmwrite([filename, 'V.txt'],[x',Va']);
    dlmwrite([filename, 'R.txt'],[x',Ra']);
    dlmwrite([filename, 'Err.txt'], Err);
end

%% Plotting data
L = [1000,2000,3000,4000,5000];
N = [2001,4001,6001,8001,10001];
linespec = {'-b', '-c', '-m','-y', '-g'};
figure; hold on;
%ax1 = subplot(2,1,1); xlim([0,200]); hold on;
ax2 = subplot(1,1,1); xlim([1,500]); hold on;
for i = 1:5
    filename = sprintf('MoS2(f=%1.1E,L=%i,N=%i)',1E-4,L(i),N(i));
    V = dlmread([filename, 'V.txt']);
    R = dlmread([filename, 'R.txt']);
    %plot(ax1,V(:,1),V(:,2), linespec{i}); 
    
    plot(ax2,R(:,1),1./R(:,2), linespec{i});
    
end
x = linspace(-1500,1500,6001);
R_metal = 2*(0.0318)/(2*pi^2)./x;
plot(ax2,x,1./R_metal,'-k', 'linewidth',2);
    
    
    
    
    